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ABSTRACT 

Theoretical models predict that a population of Intermediate Mass Black Holes (IMBHs) of 
mass 10^“^ Mq might form at high (z > 10) redshift by different processes. Such 

objects would represent the seeds out of which z ^ 6 Super-Massive Black Holes (SMBHs) 
grow. We numerically investigate the radiation-hydrodynamic evolution governing the growth 
of such seeds via accretion of primordial gas within their parent dark matter halo of virial tem¬ 
perature Tyir ^ 10^ K. We find that the accretion onto a Direct Collapse Black Hole (DCBH) 
of initial mass Mq = 10^ Mq occurs at an average rate M, 1.35 MEdd — 0.1 Mq yr“^, 
is intermittent (duty-cycle ^ 50%) and lasts 142 Myr; the system emits on average at 
super-Eddington luminosities, progressively becoming more luminous as the density of the 
inner mass shells, directly feeding the central object, increases. Finally, when ^ 90% of the 
gas mass has been accreted (in spite of an average super-Eddington emission) onto the black 
hole, whose final mass is ^ 7 x 10^ Mq, the remaining gas is ejected from the halo due to a 
powerful radiation burst releasing a peak luminosity Lp^ak ^ 3 x 10^^ ergs“^. The IMBH 
is Compton-thick during most of the evolution, reaching a column density Nh ^ 10^^ cm“^ 
in the late stages of the simulation. We briefly discuss the observational implications of the 
model. 

Key words: accretion - black hole physics - hydrodynamics - radiative transfer - dark ages, 
reionization, first stars - methods: numerical - cosmology: early Universe 


1 INTRODUCTION 


In the standard ACDM cosmological model, the formation of the 
first stars occurred in the redshift range 30 > z > 20 (|Barkana| 
& Loeb||200l| |Miralda-Escudg||200^ [Bromm & Yoshida||2011[ 


Stacy & Bromm|2014[|Hummel et al.|2014| for a recent and exten- 
sive review see also |Bromm|2013| l in mini-halos, i.e. dark matter 
structures with virial temperatures Tvir ^ 10^ K and total masses 
Mh ^10® Mq. The concurrent formation of the first black holes 
([Bellovary et al.|2()TT| Volonteri & Bellovary|20T^|Agarwal et al.| 


|2013||2014| see |Hainiian||2013| for an updated review), as a final 


product of the evolution of the first massive stars, represents a sec¬ 
ond key event during the same cosmic epoch. The appearance of 
these classes of objects likely had a very strong impact on both the 
interstellar and the intergalactic medium, due to their radiative and 


mechanical feedback (|Park & Ricotti|2011[ 

Petri et al.|2012||Park 

& Ricotti|2012[|Jeon et al.|2012[|Tanaka et a 

.|2012[|Maiolino et al. 

20121|Park & Ricotti|20131|Nakauchi et al.|2014>. 


The process of cooling and fragmenting the gas plays a role 
of key importance in this cosmic epoch. Mini-halos primarily cool 
their metal-free gas through molecular hydrogen line emission. Un¬ 
der intense irradiation in the Lyman-Werner band (LW, = 
11.2 — 13.6 eV), H 2 is photo-dissociated via the two-step Solomon 
process (see the original study by |Draine & Bertoldi|1996| l, so that 
cooling (i.e. the formation of stars and stellar mass black holes) 


is quenched ( |Visbal et al.|2014| ). The UV specific intensity in the 
LW band, jJ — J 21 x 10“^^ erg s“^cm“^Hz“^ sr“^, required 
for quenching is somewhat uncertain, but lies in the range J 21 = 
0.1 — 1 (see, e.g., |Machacek et aH2001| and |Fialkov et al.|2013] ). 
When instead primordial, atomic-cooling halos (Tvir > 10^ K) are 


exposed to a LW flux of even higher intensity, > J* (Loeb| 

& Rasio||1994| Eisenstein & Loeb||1995| |Begelman et al.||2006| 

Lodato & Natarajan|2006||Regan & Haehnelt| 

120091 IShang et al.l 

201 0(1 Johnson et al.|2012|lAgarwal et al.|2012l 

[ Latif et al.|2013a|l 


the destruction of H 2 molecules allows a rapid, isothermal col¬ 
lapse, initially driven by H i Lya line cooling, later replaced by 
two-photon emission. The precise value of J* depends on several 
factors, but there is a general consensus that it should fall in the 
range 30 < J'l < 1000, depending on the spectrum of the sources 
( |Sugimura et al.|2014] l. 


Several theoretical works (|Bromm & Loeb||2003 

Begelman 

et al. 

2006 II Volonteri et al.|2008| Shang et al.|2010| Jo 

inson et al. 

2012 

(have shown that the result of this collapse is the formation of 


a Direct Collapse Black Hole (DCBH) of mass M, 10^“^ Mq, 
likely passing through the intermediate super-massive protostar 
phase, either directly collapsing into a compact object due to gen¬ 


eral relativistic instability ( |Shibata & Shapiro|2002[|Montero et ^ 
|2012| ) or at the end of a very rapid evolution if the super-massive 
star reaches the Zero-Age Main Sequence (ZAMS, for a more de- 
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tailed discussion see [Ferrara et al.||2014] >. Once formed, the sub¬ 
sequent accretion of the remaining gas from the parent halo leads 
to a further growth of the DCBH into a fully-fledged Intermediate 
Mass Black Hole (IMBH) of mass M, 10^“® M©. This sce¬ 
nario is confirmed by cosmological simulations, as those presented 
by jLatif et aL| ( |2013^ , who have shown that under the previous con¬ 
ditions (atomic-cooling halos irradiated by Ju > J*) very strong 
accretion flows up to 1 M© yr“^ may take place. As calcula¬ 
tions by jHosokawa et al.| ( |2013| ) and jSchleicher et al^ ( |2013| ) suggest 
that the effects of radiative feedback from the accreting protostar is 
weak, due to its cool 6000 K) photosphere, it has been possible 
to safe ly extend the simulations to longer time scales ( jLatif et a^ 
2013b| ), up to the formation ofa~ lO^M© central condensatioijj 


In our work, we follow the accretion history onto a newly 
formed DCBH, with special emphasis on the dynamical and radia¬ 
tive properties of the inner parts of the parent halo, providing most 
of the accretion material to the central black hole. In particular, we 
aim at clarifying a number of aspects, including (i) the accretion 
time scale; (ii) the time-dependent emitted luminosity; (iii) the final 
outcome of the accretion phase, with a progressive depletion of the 
gas or a final outburst; (iv) the fraction of the halo gas accreted by 
the black hole and that ejected in the outskirts by radiative feedback. 

Previous works have already attempted to describe the accre¬ 
tion process onto black holes of different sizes. jSakashita & Yoko-j 
jsawaj ( |1974| ) employed the method of similarity solutions to de¬ 
scribe the time evolution of the accretion process in the optically 
thin regime (i.e. r < 1, where r is the optical depth). [Tamazawa 
jet al.| ( p^75] ), instead, used a full radiative transfer approach to de- 
scribe the process in the optically thick regime (i.e. t > 1) assum¬ 
ing steady state, i.e. without any explicit time dependence in the 
accretion flow. In contrast with these previous works, we aim at a 
full time-dependent description in the optically thick regime, sim¬ 
ilarly to the somewhat idealized approach used by ( [Park & Ricottij 
|2011||201^|2013[pohnson et al.|2013| l, but extending these studies 
in several aspects and including a more complete physical descrip¬ 
tion of the accretion process. Recent works (Alexander & Natarajanj 


|2014[[%lonteri & Silk|2014}|Madau et al.|2014| | have proposed the 

occurrence of brief, recurring but strongly super-critical accretion 
episodes (with rates even 50 — 100 times larger than the Eddington 
limit) to explain the rapid black hole mass build-up at high redshifts. 
An early phase of stable super-critical quasi-spherical accretion in 
the BHs was also proposed in jVolonteri & Ree^ ( |2005| ). Such large 
accretion rates may be sustainable in the so-called “slim disk” so- 
lution (jBegelman & Meier|1982[|Paczynski & Abramowicz|1982[ 


[Mineshige et al.|2000[|Sadowski|2009 [201 Ij l, an energy-advecting, 
optically thick accretion flow that generalize the standard thin disk 
model ( jShakura & Sunya^|1976| ). In these radiatively inefficient 
models the radiation is trapped into the high-density gas and is 
advected inward by the accretion flow: as a consequence, the lu¬ 
minosity is only mildly dependent on the accretion rate and very 
large flows are sustainable. These works, while investigating the 
accretion flow at much smaller spatial scales, offer an interesting 
perspective in the discussion of the implications of our efforts, as 
detailed in Sec.|4] 

The present work is the first necessary step towards a precise 
prediction of the observable properties of IMBHs, whose existence 
has remained so far in the realm of theoretical, albeit physically 
plausible, speculations. This effort seems particularly timely given 


^ Even higher masses can result if magnetic fields, suppressing fragmenta¬ 
tion, are included jLatif et al.|2014| . 


the advent of powerful telescopes as the James Webb Space Tele¬ 
scope (JWST), becoming available in the next few years, and high-z 
ultra-deep surveys like the HST Frontier Fields ( jCoe et al.|2014| t. 
In practice, we aim at determining the Spectral Energy Distribution 
of IMBHs in order to build diagnostic tools able to uniquely iden¬ 
tify these sources among the other high-z ones. If successful, this 
strategy would not only represent a breakthrough in the study of the 
first luminous objects in the Universe, but it would also shed some 
light on the puzzles provided by the interpretation of the formation 
of SMBHs (see e.g. jFan et al.|2001[[Mortlock et al.|201 l[|Petri et al.j 
|2012j l and the excess recently found in the power spectrum of the 
cosmic Near Infrared Background fluctuations (for an overview, see 
|Yue et al.|20I3| ). These issues will be discussed in more details in 
Sec. in 

The outline of this paper is as follows. In §2 we describe the 
physics and equations of the radiation-hydrodynamic problem we 
aim at solving, along with the numerical implementation and ini¬ 
tial conditions. In §3 we present the results of our simulations, 
providing a full picture of the accretion and feedback processes. 
Finally, in §4 we provide some discussion and the conclusions. 
The Appendix contains some more technical aspects of the sim¬ 
ulations. Throughout, we adopt recent Planck cosmological pa¬ 
rameters (|ManckConaboration[20T3]): (Qrn, CTs) = 

(0.32,0.68,0.05,0.67, 0.96,0.83). 


2 PHYSICAL AND NUMERICAL IMPLEMENTATION 


The present study is based on a series of radiation-hydrodynamic 
simulations. Our code is designed to execute a fully consistent 
treatment of uni-dimensional spherically-symmetric hydrodynamic 
(HD) equations and a simplified version of Radiative Transfer (RT) 
equations. 

The simulated spatial region is large enough to allow us to ne¬ 
glect deviations from spherical symmetry, e.g. the presence of an 
accretion disk which may form at much smaller scales. As detailed 
in jde Souza et al.| ( [2013| ), the primordial halos which generated the 
DCBHs rotate very slowly, with a mean value of the spin parameter 
A = Jc/(GM^) = 0.0184, where J is the angular momentum of 
the halo with mass Mh and c is the speed of light. Under these con¬ 
ditions, deviations from the spherical symmetry become important 
at the centrifugal radius Rc — ~ 10“® pc ~ 

100 i^s, with denoting the Schwarzschild radius: this value is 
~ 10^ times smaller than the internal boundary of our simulations. 
Our resolution is designed to resolve the Bondi radius (or gravita¬ 
tional radius, see |Bondi|1952| l: 

/is = = 1.5pc = 5 X lO'^flvir = 10®fls (1) 

^s(oo) 

where Rr^^ir is the virial radius of the halo and Cs^oo) is the sound 
speed at large distances from the accretion boundary, defined as: 


Cs(oo) 



15kms ^ 


( 2 ) 


Here, 7 = 5/3 is the ratio of specific heats, p{Rb) and p{Rb) 
are the gas pressure and mass density at the Bondi radius, respec¬ 
tively. The Bondi radius largely varies during our simulations, but, 
for clarity reasons, all distances in the plots are expressed in terms 
of its initial value RB{t = to). Interestingly, even Adaptive Mesh 
Refinement cosmological simulations cannot resolve this spatial ra¬ 
dius (see, e.g. jPallottini et al.|2013| where the maximum resolution 
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is ~ 5 kpc) and therefore they have to resort to some kind of sub¬ 
grid prescriptions for the black hole growth. In this context, the 
usual methodology to deal with black holes is to suppose that they 
irradiate with luminosities L ^ LEdd where: 


Lsdd = 3.2 X 10^ 



(3) 


is the Eddington luminosity. The domain of our simulations spans 
approximately from 0.1 Rb to 2Rb. This spatial dynamic range 
is designed to encompass with the highest possible resolution the 
spatial region of interest for the full simulation, i.e. from the radius 
of gravitational influence Rb) down to the smallest radius 
0.2 Rb) reached by the propagating density wave (see Sec.|^. This 
spatial range has been successfully tested to verify the convergence 
of the main quantities derived in this work. The natural time scale 
of the problem is given by the free-fall time at the Bondi radius: 


iff 


1 

^GpiRB) 


10" yr; 


(4) 


Here p{Rb) ~ 3 X 10 g cm ^ is the mass density at the Bondi 
radius at the beginning of the simulations. The time scale for the full 
RT simulation is about ~ 10^ t//, due to the presence of radiation 
pressure which slows down the collapse. 

In the following subsections we describe the physics included 
in the simulations, separating the HD and the RT parts for clarity 
reasons. Some more technical aspects (boundary conditions, two- 
stream approximation, heating and cooling terms and photon diffu¬ 
sion) are deferred to the Appendix. 


2.1 Hydrodynamics 

We solve the standard system of ideal, non-relativistic Euler’s equa¬ 
tions, i.e. neglecting viscosity, thermal conduction and magnetic 
fields, for a primordial (H-He) composition gas, spherically accret¬ 
ing onto a central DCBH, supposed at rest; the angular momentum 
of the gas with respect to the central object is zero. 

The code evolves in time the following system of conserva¬ 
tion laws for mass, momentum and energy, solving for the radial 
component: 

g+V-(pv) = 0 (5) 

^ + V • (q X V) = -(7 - 1)VS - Vprad +gp ( 6 ) 

^ + V • (e;v) = -(7 - 1)EV • V + (if - C) (7) 

where p is the gas mass density, v is the gas velocity (taken to 
be positive in the outward direction), q = pv is the momentum 
density and E is the energy density. Moreover is the additional 
radiation pressure, g(r) is the gravitational field generated by the 
central black hole and {H—C) are the heating and cooling terms. In 
the following, we neglect the vector notation since we only consider 
the radial components of the previous quantities. The total energy 
density E is given by the relation: 

E = pt+^ (8) 

where e is the specific gas thermal energy: 

t=-l-lRT ( 9 ) 

7 - 1 M 
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Here, T is the gas temperature, R is the gas constant and p is the 
mean molecular weight which, for a primordial gas with helium 
fraction Yp = 0.2477 ( [Peimbert et al.|2007] ) and no metals Zp = 0, 
is equal to p = 1.15. The gas thermal pressure is given by the usual 
equation for ideal gases: 

Pg = IpRT, (10) 

p 

while the gravitational acceleration is 


9{r, t) 


GM.jt) 


( 11 ) 


The value of the black hole mass M, (t) changes with time, due 
to the accretion, with the following set of rules, where M, = 
47rr^p|u|: 


rM.(t)/0 v{r = TminR) < 0 

\M.{t) = Mo + il-v)f^M.dt 

where Mo = lO" M© is the initial value for the mass we adopt 
for the DCBH, M, denotes the time derivative of M, and p is the 
efficiency factor for mass-energy conversion (see the RT subsection 
below for more details). 

The system E qs.[^- [7] are solved with a Linearized Roe Rie- 
mann solver ( |Roe||19^ , a method based on Roe’s matrix char¬ 
acteristic decomposition, which offers superior quality in treating 
the diffusion in hydrodynamic quantities. The time-stepping algo¬ 
rithm is a classical Runge-Kutta third-order scheme. The time step 
is computed from the well-known CEL condition: 

V^=G< Crnax (13) 

Ar 

where C is the dimensionless Courant number and Cmax = 0.8. 


2.2 Radiative Transfer 


In our simulations, all radiation-related quantities are integrated 
over frequencies: this treatment serves as a good approximation for 
the main radiative features. Our RT modeling builds upon the works 
by |Tamazawa et al.| ( |1975| ) for what concerns the general theoretical 
framework; we also exploit the computationally effective scheme 
by |Novak et al.| ( [2012| ). In the following we first introduce the rele¬ 
vant RT equations and then discuss their numerical implementation. 

In the usual notation, J is the intensity of the radiation field, 
S is the source function, H (K) is the first (second) moment of 
intensity. All these quantities are functions of time and position. 
The closure relation between the second moment of intensity K 
and the intensity J is given by the so-called Eddington factor T 
( [Hummer & Rybicki|1971[[Tamazawa et al.|1975| l, here defined as: 


K ^ 1 + t/tq 
J 1 + 3t/to 


(14) 


where tq is the optical depth at which the Eddington factor becomes 
equal to 1 / 2 . 

The source total luminosity L and H are related at any point 
of the spatial domain by L = In turn, L depends on the 

gas accretion rate M, onto the IMBH (for u < 0): 


L, = pc^Mm = pc^ (47rr^p|u|) (15) 

where we employ an efficiency factor p = 0.1 (see e.g. |W&| 
|Tremaine|2002[ [Madau et al.|2014^ . Generally speaking, the effi- 
ciency factor ranges from p = 0.057 for a Schwarzschild (i.e. non¬ 
rotating) black hole to 77 = 0.32 for a maximally rotating object 
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(see |Thorne|1974^ . We then define L» the luminosity evaluated at 
the innermost grid cell, located at r = Tmin'- 

j |'u(r7TT,in) Ij if 0 /i^\ 

= 0 ifv(rmin)>0 

In addition, we define fEdd = M^jMEdd'- the accretion is super¬ 
critical (super-Eddington) if fEdd > 1- Note that, only in the case 
of a fixed value of the efficiency factor 77 , fEdd = Mm/MEdd = 
LjLEdd holds. The acceleration caused by radiation pressure is 
then: 


_ n{p,T)L 

^rad — AO 

47rr^c 


(17) 


where A^(p, T) is the total opacity of the gas: in our simula¬ 
tions the Thomson electron scattering, bound-free and H~ opac¬ 
ity terms are included. The Thomson opacity is calculated follow¬ 
ing the temperature-dependent prescription given in |Begelman et al.| 
( |2008| ), namely: 

kt { t ) = 0 . 2(1 + m ) ^ ^ ( r ) r ,)-/3 g ~'’ 

with P = 13. Below T ^ — 8000 K, the opacity rapidly falls 

due to the decrease of the ionized fraction: as a consequence, also 
the effectiveness of the radiation pressure is quenched. The bound- 
free opacity is computed through the Kramers’ approximation for a 
metal-free gas (Z = 0) (see e.g. [Carroll & Ostlie|1996| ): 

K6/(/5,T) = 4 X 10""(1 + Xp)/5T-^'"(1 - Xion) (19) 


Car- 


where Xp is the primordial hydrogen fraction. The opacity d ue to 
the negative hydrogen ion H~ is computed as (see again e.g. 

[roll & Ostlie|1996| ): 

Hh- {p, T) = 3.5 X 10-"- Xion) (20) 


With a metal-free gas {Z — 0), within the ranges of mass densities, 
temperatures and ionized fractions in our spatial domain, the dom¬ 
inant source of opacity is the Thomson one. Our definition of the 
optical depth r is: 

t{R) = — f Kpdr (21) 

Jo 

The quantity reported in the graphs is the optical depth for an exter¬ 
nal observer. 

As we will see, due to feedback effects, accretion onto the cen¬ 
tral object occurs in an intermittent manner. It is then useful to intro¬ 
duce the duty-cycle, defined as the fraction of time spent accreting 
within a given time frame of duration Ttot - 

V = XiSSL ( 22 ) 

Ttot 


The instantaneous value of the duty-cycle is computed dividing the 
total integration time Ttot in slices and computing the duty-cycle 
with respect to each slice. 

The equations for steady and spherically-symmetric transfer of 
radiation have been derived, e.g. by [Castor] { |1972| ) and [Cassinelli &[ 
[Castor](T973[ ). For our problem, it is appropriate to assume steady- 
state RT equations since the light-crossing time at the Bondi radius, 
Rb/c ^ 5yr, is negligible with respect to the free-fall time, see 
Eq.|^ The full RT equations are: 


c dr \ p J c dr \p J P c \ r J 


1 dL 
Airpr^ dr 


47r/^( J — S) 


(23) 


dK ^ r 3K-J \ ^ vdH 
dr \ r J c dr 


- --^H 
r c p c dr 


phvH 


( 24 ) 


The term k,{J — S) handles the gas heating and cooling, as the 
{H — C) terms in the energy equation, see Eq.[7] These equations 
are correct up to the first order in jd — v/c and are suitable for high¬ 
speed accretion flows, where fd is not negligible. In the previous 
equations, the transition from the optically thin to the optically thick 
regime is described by the density-dependent Eddington factor T 
(decreasing from 1 to 1 /3 with increasing optical depth) and by the 
term ~ p/^(p, T). 

[Novak et al.[ \2012\ presented several computationally- 
effective non-relativistic RT formulations which yield the correct 
behavior both in optically thin and optically thick regimes. These 
formulations are convenient because they allow to obtain accurate 
results with a lower computational complexity. We defer the inter¬ 
ested reader to the full derivation in [Novak et al.[ {2012\ and we 
write down only the final formulation. 

^ = A-nFiE - AirpKj) (25) 


dJ 2Jw (3 — 2w)pkL 

dr r IbTT^r^ 


(26) 


In these expressions, w is an interpolation parameter which controls 
the transition from optically thin to optically thick regimes and is 
a function of the position: it ranges from w — d (for an isotropic 
radiation field, i.e. in the optically thick regime) to = 1 in the 
optically thin regime. Furthermore, E is the source term, playing 
the same role of S in the relativistic treatment (see the Appendix 
for a m ore det ailed description). Eq. [^is derived from the full 
RT Eqs. 23|24 by neglecting the u/c terms and using E instead of 
S. Eq. [26] instead, is derived by using the specific values of the 
Eddington factor E = Kj J and the interpolation factor w in each 
regime: {E — 1, w — 1) for the optically thin and {E — 1/3, 
u; = 0 ) for the optically thick. 

Finally, we need to specify the boundary conditions. As in the 
study of stellar interiors, the second order differential equation for 
L(r), or the two first-order ODEs in L(r) and J(r) (Eqs. 
requires two boundary conditions, at the inner and outer bound¬ 
aries. The innermost cell of the grid radiates the luminosity pro¬ 
duced by the accretion flow onto the black hole (see Eq. [^, so 
that L{rmin) = T,. Far from the black hole, the luminosity is ex¬ 
pected to resemble a point source because the scattering becomes 
negligible, so that: 


J, . ^ L{rmax) ,271 

J{rmax) - 2„2 > 

In order to solve the set of boundary-value differential Eqs. \25\ 
\26\ the so-called s hooting method with Newton-Raphson iteration 
( Press et al.|1992[ ) works well up to the optical depth of a few. Be¬ 
yond this limit, the shooting method becomes unstable and it is 
necessary to switch to a more powerful relaxation method ( [Press] 
[et al.[[T9^ . However, following the evolution of the system for 
physically relevant time scales (~ 100 Myr) requires such a large 
number of steps that even the relaxation method becomes computa¬ 
tionally not viable. To overcome this problem, we follow the two- 
stream approximation method outlined in the Appendix B of [Novak[ 
[et al.[ ( | 20 T^ and sketched in our Appendix. 
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Table 1. Simulation parameters for the two test simulations (T1 and T2) and for the full one (FS) which includes a complete solution of the radiative 
transfer. 


Parameter T1 T2 FS 


'f'min [pc] 
'f'max [pc] 

Integration time [yr] 
Radiation pressure 
Energy equation 


0.16 

3.2 

3.2 X 10^ 
no 

adiabatic 


0.16 0.16 

3.2 3.2 

3.2 X lO-^ 1.4 X 10® 
yes yes 

adiabatic non-adiabatic 


2.3 Initial Conditions 


We model the spherically-symmetric gas accretion onto a seed 
black hole of initial mass Mo = 10^ M©, assumed at rest at the 
center of a dark matter halo of virial temperature Tvir ~ 10^ K 
and total mass (dark matter and baryons) Mh = 6.2 x 10^ M© at 
redshift z = 10. For r Rvir most of the mass is baryonic: for 
this reason we supposed that the gas mass contained in our compu¬ 
tational domain (r < 2Rb Rvir) is Mgas — {^b/^m) Mh = 
9.6 X 10® M© (similarly to what has been done in 


Latif et al. 


2013a 


see their Fig. 5). This assumption is reasonable since, during the 
formation process of the DCBH, most of the baryonic mass of the 
halo is expected to fall into its gravitational influence (the halo has 
no spin, then no centrifugal barrier is formed). 

We assume that the gas follows the density profile derived 
from the simulations in [Latif et al.| ( |2013^ , which is well approxi¬ 
mated by the following functional form: 


p{r) 


Po 

1 +(r/a)2 


(28) 


where a is the core radius, estimated from Fig. 1 of [Latif et al.j 
( 2013a i as a ~ 270 AU ~ 1.3 x 10“^ pc: this scale is too small to 
be resolved by our simulations. Normalization to the gas mass con¬ 
tent, Mgas, gives the central density value po ~ 10“^^ gcm“^. 
From the initial prescription for the density held and assuming an 
isothermal profile with T = 10^ K (as in Latif et ar]|2013a and 
leading to a very small value, ~ 10“^, of the ionized fraction 
throughout the domain), the initial conditions for the pressure held 
are derived directly from the equation of state, Eq. The initial 
conditions for the radial velocity are set to a very small, spatially 
constant value vq < 0. After a very brief (<C t//) transient, the 
system adjusts to a velocity profile consistent with a rapid accretion 
across the inner boundary of the grid. In addition, the initial density 
profile is also rapidly modified within the Bondi radius, while out¬ 
side the gravitational influence of the black hole it remains almost 
unaltered. 


The spatial grid for all simulations extends from Tmin — 
5.0 X 10^^ cm 0.16 pc to Vmax — 1-0 X 10^^ cm 3.2 pc, 
with 3000 logarithmically spaced bins. This range is optimized for 
our computational target, allowing to follow the entire displacement 
of the density wave and extending out to the Bondi radius. In our 
work we do not have the sufficient spatial resolution to resolve the 
HII region produced by the accretion-generated radiation (but see 
e.g. [Milosavljevic et al. [2009] where the HII region is resolved). A 
simple estimate of the extension of the Stromgren radius in our case 
yields a dimension of ~ 0.03 pc, which is in very good accordance 
with the radius reported in [Park & Ricotti[ ( [2011[ t. 

For further details about the simulations, e.g. the boundary 
conditions, see the Appendix. 


3 RESULTS 

We present our results starting from the two test simulations T1 
(pure hydro) and T2 (hydro + radiation pressure) in order to 
highlight some physical key features of the problem. Next, we turn 
to the analysis of the full simulation, FS, that in addition includes 
the complete solution of the radiative transfer; the FS simulation 
contains our main findings. The simulation parameters used in the 
three runs are reported in Table 

The total integration time Ttot for the runs TI and T2 has 
been chosen in order to show the most important features of their 
radiation-hydrodynamic evolution, while for the FS simulation it 
corresponds to the total history of the system. 


3.1 Test simulations 

3.1.1 Tl: Pure hydro 

Simulation TI is a purely adiabatic hydrodynamic simulation. The 
only two forces acting on the gas are produced by the gravitational 
held of the black hole and by the pressure gradient; therefore, it 
corresponds to the “classical” Bondi solution ( [Bondi[[1952[ ), with 
the only exception consisting in the limited gas reservoir, which 
prevents a genuine steady state to take place. The spatial profiles 
predicted by Tl are shown in Fig.[^ 

The evolution occurs over a free-fall time scale, iff ~ 3.2 x 
10^ yr during which the system progressively approaches the Bondi 
solution, reported on the density panel as a black dot-dashed line, 
in which the gas density scales as: 

j (29) 

This scaling cannot be sustained for t iff, due to the fact that 
the gas reservoir is limited (this is not the case in the classical Bondi 
solution) and ceases to be valid near the Bondi radius, where the 
gravitational influence of the black hole terminates. 

The system is progressively emptied from the inner part of the 
environment, due to the absence of radiation pressure: the density 
near the inner boundary drops by a factor ~ 5 during the simulation 
time and the effect is propagated outward, up to the Bondi radius. 
The velocity of the gas increases with time as well, stabilizing to 
a value of order ~ — 100kms“^ at the inner boundary, which is 
the result of an equilibrium between the gravitational acceleration 
and the decreasing thermal pressure of the gas. The temperature of 
the infalling gas increases, reaching peaks of 4 x 10^ K in the in¬ 
ner regions. The temperature profile is reflected in the behavior of 
the ionized fraction (see the Appendix for the relevant equations), 
which is smaller than 1 only when the temperature drops below the 
~ 2 X 10^ K level. The small jumps visible at the outer boundary 
of the temperature spatial profile are due to the different rates of 
change of density and pressure (remind that T oc p/p, as in Eq. 
P^ . The gas, outside the sphere of gravitational influence of the 
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Figure 1. Spatial profiles for the T1 simulation: total integration time is Ttot ~ 3.2 x 10^ yr. The purple line labeled as “IC” represents the initial conditions: 
where this line is not present the corresponding initial condition is impossible to show on the plot. The colored lines correspond to different times of the 
simulation, equally spaced such that Ttot = i At with At = 3200 yr and i = 1, ... 10 (for clarity reasons, only the labels i = 1 and z = 10 are shown on 
the plot). All horizontal axes are in units of Rb- The density panel reports the classical Bondi solution with a black dot-dashed line. The accretion rate is the 
mass flux at a given radius and it is plotted with the velocity sign, i.e. a positive (negative) value corresponds to an outflow (inflow). 


black hole, is swept away by the thermal pressure (dp/dr < 0 in 
the entire spatial domain, so that the related force is always in the 
outward direction) causing an abrupt change in the spatial scaling 
of pressure and density. This effect is reflected most dramatically on 
the ionized fraction, since this quantity is very sensitive to modifi¬ 
cations of the temperature around the level ~ 10^ K. These effects 
are unimportant for the overall evolution of the system and are vis¬ 
ible also in the T2 simulation. 

The flow accretes at strongly super-Eddington rates at all 


times, reaching peaks of fsdd = M*/Msdd ~ 2000: the Edding¬ 
ton limit does not apply in this case due to the absence of radiation 
pressure. The accretion rate progressively decreases due to the den¬ 
sity drop: this causes the slow reduction of the luminosity at the 
innermost cell, clearly visible in the bottom-right panel of Eig.[2 
The emitted luminosity is obscured by high column densities 
at the beginning of the simulation, while the drop of the density 
decreases the optical depth with time, down to a value ~ 6. The 
spatial profile of the emitted luminosity is described by Eq.[^(and, 
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more specifically, by Eq.j^in the Appendix). The IMBH mass at 
the end of the T1 simulation is M, 2.0 x 10^ M©. 

3.1.2 T2: Adding radiation pressure 

In the T2 simulation an outward radiative force, corresponding to 
a fixed (i.e. not tied to M.) value of the luminosity L = 2 x 
lO^^ergs”^, is added to the gravitational force and to the pres¬ 
sure gradient, while the energy equation is still adiabatic. The ra¬ 
diative force is active only when the black hole is accreting, i.e. 
when v{rmin) < 0. The aim of the T2 simulation is to show in a 
simplified way the effect of the radiation pressure on the gas. 

Fig .[^is a comparison between the fixed value of the luminos¬ 
ity employed in T2 and the luminosity resulting from the T1 simu¬ 
lation, computed through the accretion rate M, at r = Tmin with 
Eq.[^ but not included in T1. The latter is ~ 3 orders of magnitude 
larger, due to the absence of radiation pressure quenching the accre¬ 
tion flow. The Eddington luminosity is also shown for comparison, 
its progressive increase being due to the change in M,. The value 
of L is set in order to be at any time larger than the corresponding 
LEdd- 

The radiation pressure can modify the accretion flow in two 
ways. If fEdd = MmjMEdd < 1, the effect is a decrease of the 
accretion rate M,. Naming Mti the accretion rate without any ra¬ 
diative force (i.e. the one in the TI simulation) and Mr the accre¬ 
tion rate with the addition of a sub-Eddington radiation pressure, it 
is easy to show that the following relation holds: 

Mr = Mti V^l — fEdd (30) 

If, instead, fEdd > 1, the accretion is intermittent (i.e. V < 1): 
the infalling gas is swept away by the radiation pressure and some 
time is needed to re-establish the accretion. Under the simplifying 
assumption that the initial infalling velocity of the gas is slow, it is 
possible to show that, if fEdd > 1, an estimate of the value of the 
duty-cycle is given by: 

■D = {2fEcid - 1)"^ (31) 

Under these assumptions, P = 1 for fEdd = 1, while for fEdd > 1 
it steadily decreases. 

The T2 simulation is an example of the latter case, with 
fEdd ~ 1-5. From the previous analysis we expect two major dif¬ 
ferences with respect to the Tl simulation, namely: (i) the IMBH 
accretes ~ 50% less mass (if the total integration times are equal) 
because V is smaller by a factor ~ 2, and (ii) the IMBH pro¬ 
duces some feedback effect on the surrounding gas. This is ex¬ 
actly what we observe in the T2 simulation. The final mass is 
M, = 1.4 X 10^ M©, i.e. the black hole accreted ~ 60% less 
mass than in the Tl simulation, in agreement with our rough esti¬ 
mate in Eq. [ST] In addition, the spatial profiles for this simulation, 
shown in Fig.l^ manifestly provide some evidence of the effect that 
the radiation pressure exerts on the gas. A density wave, produced 
by the radiative feedback, propagates in the outward direction with 
velocities up to ~ 20 km . It is interesting to note that this wave 
is mildly supersonic, since the value of the sound speed in a gas at 
T ~ 2.6 X 10^ K is ~ 19kms“^. In addition, the positive val¬ 
ues of the accretion flow measured in a large part of the spatial 
domain indicate the occurrence of a gas outflow from the exter¬ 
nal boundary. The high-density wave, with temperatures as high as 
~ 2.8 X 10^ K, is followed by a rarefaction zone, where the temper¬ 
ature drops to ~ 7000 K (the temperature profile follows the pres¬ 
sure, i.e. the cooling is adiabatic), decreasing the ionized fraction to 
very small (~ 10“^) values as well. The accretion flow promptly 



Time [yr] 

Figure 2. Comparison between the emitted luminosity for different simula¬ 
tions (note that the vertical axis is broken). The red line is the luminosity of 
the Tl simulation (computed from M* at r = rmin with Eq.[^ but not 
included in the physics). In green, the fixed luminosity value used for the 
T2 simulation L = 2x 10^^ erg . In blue, for comparison, the Edding¬ 
ton luminosity for the T2 simulation, whose value increases along with M*. 
The value of L is always above the Eddington limit. 


(after ~ 6000 yr) stabilizes to a value fEdd ~ —500 at the in¬ 
nermost cell: this value is set by the fixed radiation pressure of the 
T2 simulation, which indirectly sets also the velocity at which the 
accretion flow is re-established at the end of each idle phase (the 
larger is the radiation pressure, the longer is the time needed for 
accretion to re-start, then the larger is the resulting mass inflow). 

This general framework is explained by the following sce¬ 
nario: when the black hole accretes, the fixed super-critical emitted 
luminosity sweeps away the surrounding gas, affecting a spherical 
region of radius rr, defined such that T(rr) = 1. The gas located at 
r ^ Tr is also accelerated upward, not by the radiation pressure in 
this case but by the thermal pressure, and acquires a positive veloc¬ 
ity. When the irradiation is temporarily shut down, the gas located 
at r Tr is affected by the strong gravitational field of the black 
hole and falls back in due course. 

3.2 The Full Simulation 

The aim of the FS simulation is to describe the accretion flow onto 
a DCBH with initial mass Mo = 10^ M©. The final integration 
time for this simulation, when all the gas contained in the halo is 
expelled by radiation pressure, is ~ 142 Myr. The forces acting on 
the gas are the gravity of the black hole, the pressure gradient and 
the radiation pressure. 

The differences with respect to the previous Tl and T2 simu¬ 
lations are: (i) the radiation pressure is computed self-consistently, 
i.e. with Eq.[^and (ii) the energy equation is non-adiabatic, with 
the inclusion of bremsstrahlung cooling and atomic cooling (see the 
Appendix). 

3.2.1 Spatial structure and time evolution 

Broadly speaking, the simulation allows the identification of three 
distinct evolutionary phases of the system, described in turn below. 

(i) Initial Transient Phase: the gas, initially almost at rest (as 
detailed in Sec.[^, is accelerated downward, progressively increas¬ 
ing M, and, as a consequence, the emitted luminosity L,, as shown 
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Figure 3. Spatial profiles for the T2 simulation: total integration time is Ttot ~ 3.2 x 10^ yr. A fixed radiative force, determined by the luminosity L = 
2 X 10"^^ erg and active only when v{rmin) < 0? is added to the gravity and to the pressure gradient, as explained in the text. The colored lines correspond 
to different times of the simulation, equally spaced such that Ttot = i with At = 3200 yr and z = 1, ... 10 (for clarity reasons, only the labels i = 1 and 
z = 10 are shown on the plot). 


in Fig. 1^ This plot shows that the emitted luminosity increases by 
~ 3 orders of magnitude in only ~ 200 yr, a fraction ~ 10of 
the full evolution of the system. This process is self-regulated, due 
to the interconnection between gravity, accretion rate and radiation 
pressure: the gravity tends to increase the accretion rate by acceler¬ 
ating the gas downward, while the emitted luminosity acts against 
the infall by providing an outward acceleration. This evolutionary 
phase lasts until the emitted luminosity becomes comparable to the 
Eddington limit (see Eq.[^, approximately 10^^ erg at the be¬ 
ginning of the simulation. Above this threshold, the radiation pres¬ 


sure is able to sweep the gas away from the inner boundary and the 
accretion process becomes intermittent (V < 1). 

This initial phase lasts ~ 200 — 300 yr: the emitted luminosity 
and the fractional mass accreted (AMt/Mo — [Mt — Mo)/Mo) 
are shown in Eig|^ An estimate of the duration Tt of this initial 
phase is easily computed from the following argument. We request 
that during Tt the luminosity L, = r]c^Mm becomes equal to the 
Eddington luminosity: 


2 dM _ AirGrripC 

L^VC -TT = ^Edd = - —M 

at ctt 


(32) 
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By means of integrating between t = 0 when M(t) = Mo and Tt 
when M(t) — Mt and calling AMt = Mt — Mq we obtain: 


r]caT A AMt \ ^ ijcar AMt 

AirCmp \ Mo / AirCmp Mq 


(33) 


where the last approximation is valid for AMt/Mo ^ 1. This 
equation, with the value AM/Mq ^ 7 x 10“® taken from Fig. 
1^ gives the expected time scale Tt ~ 300 yr. Defining tEdd = 
{car)/{^T^Grrip) the Eddington time scale, the previous formula 
becomes: 


' TfiEdd 


AM 

'M^ 


(34) 


Interestingly, if we suppose that M(t) oc Mo (see e.g. 
|Silk|2014l[Madau et al.|2014l l, as in: 

'l-T] 


Volonteri & 


M(t) = Mo exp 


T] 


t 

tEdd 


(35) 


the time scale Tt is strictly independent on the initial black hole 
mass Mo. 

It is relevant to investigate in detail the mechanism which leads 
the system from a stable accretion at L ~ LEdd to the unstable 
and intermittent phase shown in Fig. for t > 270 yr. When the 
black hole starts to accrete at M, ~ MEdd, the gas near the in¬ 
ner boundary is swept away, so that the accretion is temporarily 
interrupted and the radiation pressure is turned off. The emitted lu¬ 
minosity does not affect the outer parts of the domain, for r ^ rr, 
see the similar discussion for the T2 simulation. The gas located 
in this section continues its infall, counteracted only by the ther¬ 
mal pressure exerted from the internal layers, and eventually feeds 
the black hole with increasingly larger mass infiows. The difference 
with respect to the T2 simulation is due to the direct dependence of 
the emitted luminosity from the mass infiow in the FS case. More 
specifically, this mechanism is schematically explained in Fig. 01n 
the first panel, the black hole is accreting mass from the innermost 
shell, which is collapsing just as the outer one. When the black hole 
irradiates with L > LEdd, the radiation pressure acts only on the 
innermost shell (supposing for simplicity that the outer shell is lo¬ 
cated at r ^ Tt) which is swept outward, while the outer shell con¬ 
tinues the infall. During this period, the black hole is not irradiating. 
The innermost shell eventually terminates its outward movement, 
due to the gravitational pull of the black hole. The innermost shell 
merges with the outer one, so that its overall density is increased. 
Eventually, the merged shells approach the accretion boundary and 
re-start the cycle with a larger accretion fiow, i.e. with the emission 
of a higher luminosity. 

(ii) Main Accretion Phase: this phase lasts ~ 142 Myr and is 
characterized by a duty-cycle V ~ 0.48 and an average accretion 
rate fEdd — 1.35: the accretion is super-critical on average. This 
value is computed as a global average of fEdd over the entire inte¬ 
gration time, including the idle phases (when the accretion does not 
take place) and it is in substantial agreement with the approximated 
(i.e. valid for small infiow velocities) relation given in Eq.[^ 

Fig.j^shows the spatial profiles for the FS simulation. From their 
analysis, three main features are evident: (i) a density wave ap¬ 
proaching the center, (ii) strong waves moving towards larger radii, 
visible from the velocity profile and (iii) the progressive emptying 
of the outer regions (discussed below). 

The density wave, with over-densities as high as ~ 1 order of 
magnitude with respect to the surrounding volume, moves in the 
inward direction, contrarily to the T2 simulation. This difference 
is due to the fact that in the FS simulation the radiation pressure 
and the accretion rate are interconnected, while in the T2 run the 



Figure 4. Luminosity emitted and fractional mass accreted (AM^/Mq = 
{Mt — Mq)/Mq) during the initial phase, lasting ~ 200 — 300 yr. The 
corresponding Eddington luminosity is also shown, for comparison. 



Figure 5. Schematic description of the mechanism which progressively in¬ 
creases the emitted luminosity. Here we consider only two mass shells, but 
the system must be thought to be composed of an infinite number of them. 
The black hole is at the center of each of the four panels, in black. Brown cir¬ 
cles are mass shells: the darker the shell, the higher the density. The smallest 
circle is the accretion boundary: it becomes orange when the black hole ir¬ 
radiates. Green arrows indicate accretion through the inner boundary, red 
arrows indicate a movement of the mass shell. The black dot-dashed line 
indicates the radius rr. See the text for a detailed description of the panels. 


former is fixed: this joint evolution leads to a smooth increase in 
the emitted luminosity and to a different response of the system. 
When the IMBH accretes at super-critical rates, the emitted radia¬ 
tion promptly interrupts the mass infiow and consequently the radi¬ 
ation pressure. The accelerated gas moves in the outward direction 
creating a pressure jump of a factor ~ 5 between the two sides of 
the discontinuity front (see the pressure spatial profile in Fig. [^. 
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Figure 6. Spatial profiles for the FS simulation: total integration time is Ttot ~ 1.4 x 10® yr and Ttot = i with At = 35 Myr and z = 1, ... 4 (only the 
labels i = 1 and z = 4 are shown). The luminosity lines are present only at data dumps when the IMBH is emitting. qOOO RAS MNRAS 000 000 000 
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Eventually, the gravitational acceleration of the IMBH inverts the 
velocity and the accretion starts again. 

The radiation pressure affects only a small volume, in the inner 
section of the accretion flow where r ~ 10 — 100, as the opti¬ 
cal depth spatial profile shows: the internal layers (r <C Vr) of 
the gas distribution are intermittently reached by the radiation pres¬ 
sure, while the external layers (r rr) are in a quasi free-fall 
state. For this reason, similarly to the mechanism detailed in Fig.[^ 
the density wave progressively moves towards the center, leading 
to an increase of ~ 1 order of magnitude in the density measured 
at the innermost cell, as Fig. [7] shows. The top of the density wave 
Rdw moves inward with the time scaling: Rdw oc The op¬ 

tical depth of the inner regions is increased along with the density: 
this additional effect progressively decreases the volume where the 
radiation pressure is effective. The luminosity panel is described in 
a separate subsection below. 

The velocity spatial profile shows that the outer regions (R > 
0.5Rb) are swept by waves of high-speed (10 — 20kms“^) gas. 
Note that in some regions of the spatial domain, where the tem¬ 
perature drops near ~ 10^ K, weak shocks are produced in the 
flow, with Mach numbers M ~ 1.3 — 2.0. This volume, while 
not affected by radiation pressure, is strongly affected by the ther¬ 
mal pressure exerted from the internal layers: the pressure spatial 
profile shows, in the external regions, pressure jumps as high as 
~ 6 — 7 orders of magnitude. As the radiation pressure is inter¬ 
mittent, the net result is formation of waves in the surrounding gas, 
whose magnitude and frequency increase with time: an always in¬ 
creasing energy is transported outward with this mechanism. 

Finally, the temperature spatial profile shows values as high as 
~ 10^ K in the proximity of the inner boundary, at late stages of 
accretion, caused by the very large pressure. The temperature spa¬ 
tial profile is refiected in the ionized fraction: the ionized volume 
expands outward with time. 

The complicated behavior of the mass fiux spatial profile is a 
symptom of chaotic motions occurring in the environment, caused 
mainly by the intermittent irradiation. The last data dump of the FS 
simulation shows a very small mass flux in the external layers, due 
to the fact that M, (r) (x p(r). 

(iii) Final Burst: this events marks the end of the accretion flow 
onto the IMBH, ~ 142 Myr in the simulation. The gas is swept 
away and the accretion rate goes to zero: the black hole becomes a 
dark relic, having exhaled its “last gasp”. 

As the central density increases, the emitted luminosity rises as 
well, as shown in Fig. standing always a factor ~ 2 — 3 higher 
than the Eddington luminosity. This trend is highlighted by the 
shaded region in the same plot, which shows the value of fEdd'- 
after an initial transient period lasting ~ 20 Myr (caused by the 
necessity to stabilize the accretion flow), the gap with respect to 
the Eddington luminosity increases with time: the emission is pro¬ 
gressively more super-critical. The value of fEdd here is computed 
including the idle phases: being the duty-cycle ~ 0.48, the value of 
fEdd computed considering only the active phases would be dou¬ 
bled. 

The radiated luminosity reaches the value ~ 3 X 10^^ ergs ^ 
(fEdd = M.jMEdd = LjLEdd ~ 3) and the density in the ex¬ 
ternal layers drops: the remaining mass inside the computational 
domain is a factor ~ 25 lower than the initial one, due to both 
the accretion onto the compact object and the outflow (see below). 
Differently with respect to the T2 simulation, when the radiation 
pressure was fixed to a super-critical value, its interconnection with 
the mass infiow progressively voids the halo outside-in. The exter¬ 
nal layers are not able to exert a sufficient pressure to contain the 



Figure 7. Time evolution of density and velocity computed at the innermost 
cell of the spatial grid. It most clearly shows, along with the following Fig. 

the final evolutionary phase. At the time ~ 142 Myr, a jump of a factor 
~ 3 in the velocity marks the final act in the evolution of the system: the 
remaining gas is eventually ejected outward. 


expansion of the radiation-driven shell and, as a consequence, the 
remaining gas is ejected from the system. 

This effect is most clearly visible in Fig. [7] which shows the 
velocity evolution measured at the innermost cell. The general in¬ 
creasing trend (~ 10“^ kms“^ Myr“^) of the central velocity is 
abruptly changed by a jump of a factor ~ 3. The internal layer 
starts to move outward with velocities of ~ 0.5km/s. The same 
effect is hinted at in Fig. where the accretion (and consequently 
the emitted luminosity) is abruptly interrupted. After a transient pe¬ 
riod, the velocity should be re-inverted, but the very low value of 
the gas mass still inside Rb strongly suggests that the evolution 
time scale of this system with Mo = 10^ M© is indeed of order 
150 Myr. 

5.2.2 Black hole growth 

After having investigated the space and time evolution of the accre¬ 
tion fiow, we devote some further analysis to the black hole growth, 
more specifically to the accretion time scale and to the final mass 
balance. 

An important diagnostic quantity for the accretion process is 
the duty-cycle, defined in Sec.|^ A direct comparison between the 
three simulations shows that the evolutionary time scales are quite 
different. As a simple estimator. Table reports for each simula¬ 
tion the time t 2 Mo needed for the black hole to double its initial 
mass, along with the average accretion rate < M >. For instance, 
it is instructive to compare this time scale for the T1 simulation 
(~ 10^ yr) and for the FS simulation (~ 5 x 10® yr): the aver¬ 
age duty-cycles of these simulations are different, i.e. the FS sys¬ 
tem spends a smaller fraction of time accreting. The duty-cycle is 
strictly dependent on the magnitude of the radiative force, as we 
have shown in the subsection dedicated to the T2 simulation with 
simple analytic arguments. In the FS simulation the duty-cycle sta¬ 
bilizes to an average value ~ 0.48, in agreement with the prediction 
of the approximated Eq.[^(see also |Park & Ricotti|2012| l. 

Furthermore, we have calculated the quantity of gas accreted 
by the IMBH and the amount ejected from the system, before the 
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Table 2. Diagnostic quantities for the two test (T1 and T2) simulations and for the full one (FS), specifically the time f2Mo needed to double the initial 
mass of the black hole and the average accretion rate < M >. The final mass for the FS simulation is ~ 7 x 10® Mq . 


Parameter T1 T2 FS 

t2Mo [yr] ~ 3 X 10"^ ~ 4 X 10"^ ~ 5 x 10® 

<M>[M@yr-l] ~ 3.1 ~ 1.3 ~ 0.1 



Figure 8. Time evolution of the emitted luminosity and of fsdd computed 
at the innermost cell. The smoothness of the lines is due to an averaging pro¬ 
cess over a time much longer than the typical idle periods of the accretion. 
The value of fEdd is computed as a running average over a window period 
much longer than the typical duration of the duty-cycle and includes the idle 
phases (while the average values of the luminosity do not). The blue dashed 
line shows the corresponding time evolution of the Eddington luminosity, 
which increases as the black hole mass grows. At the time ~ 142 Myr the 
accretion is abruptly terminated by the final burst. 



Time [Myr] 

Figure 9. Final mass balance for the FS simulation. All masses are normal¬ 
ized to the initial value of the gas mass, Mgas ~ 9.6 x 10® M©. The 
dashed purple line is only a reference for the unitary value. The IMBH mass 
line is smoother and extended to a longer time than the others because the 
corresponding output value is saved at a very high frequency, while other 
quantities are more discretized in time. 


final burst. The density spatial profile in Fig. shows that during 
the late stages of the simulation the outer layers of the volume are 
almost empty: in about ~ 120 Myr (the time of the last complete 
data dump) the density drops by ~ 7 orders of magnitude. The 
matter is partly accreted by the central object, partly ejected from 
the outer boundary of the system by high-speed waves, as described 
above. Our final results for the mass balance are summarized in Fig. 

The baryonic mass of the halo is reduced by a factor ~ 25 from 
the beginning of the simulation. Most of this mass (~ 90%) is ac¬ 
creted onto the black hole, in spite of an average super-Eddington 
emission, while ~ 10% is ejected by outflows. Starting from a 
DCBH of mass Mq = 10^ M©, the final mass of the black hole 
is M. - 7 X 10® Mo, a fully fledged SMBH. 


3.2.3 Radiation emission 

The study of the luminosity emitted from the IMBH is one of the 
main objectives of this work. Due to the high values of the opti¬ 
cal depth (see its spatial profile in Fig.[^, the luminosity is almost 
completely obscured for an external observer. The column density 
reaches a value of 1.3 x 10^® cm“^ in the late stages of the FS sim¬ 
ulation, as can be roughly estimated from the line number 4 in the 
density panel of Fig. considering a mean value 10“^^gcm“® 
for the mass density (10^ cm“® for the number density) between 
O.IRb and 0.577^. 


The spatial profile of the emitted luminosity is shown in Fig.|^ 
for the two data dumps when the black hole is irradiating. The lu¬ 
minosity emitted via bremsstrahlung and atomic cooling is, in this 
spatial range, completely negligible. The lines show two important 
facts: (i) the luminosity emitted near the accretion boundary slowly 
increases, due to the mechanism already detailed in this Section and 
(ii) the radiation which escapes from the outer boundary decreases 
and is blocked at progressively smaller radii. The latter effect is due 
to the accumulation of matter at lower radial distances, having in 
mind that dLjdr oc —pnL, see Eq. 37 in the Appendix. In fact, the 
density measured at the innermost cell increases by a factor ~ 10 
during the time evolution of the system, as Eig. [^demonstrates. To 
conclude, the luminosity emitted during the accretion process onto 
an IMBH at high-z is obscured for most of the system’s evolution. 
We predict that it might be observable during the final burst of ra¬ 
diation. 

In order to have a rough estimate of the observability of the 
latter phase, we suppose that the peak luminosity Lpeak ~ 3 x 
10^® ergs“^ is emitted in a small (AA/A ^ 1) range of Ear-IR 
wavelengths centered at A = 2 pm with a flat spectrum (for a de¬ 
tailed study of the contribution of DCBHs sources to the cosmic 
infrared background, see |Yue et al.|2013] l. Erom the present study 
it is hard to estimate the duration of the final burst, due to the lack 
of the necessary time resolution; however, we can put a solid lower 
limit of about 5 yr. This aspect needs to be investigated more thor¬ 
oughly by future work. Eor a source located at z = 10, the radiation 
intensity would be X 10~® Jy, which is observable by the future 
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JWST with only ~ 100 s of integration, yielding a Signal-to-Noise 
raticEJof ~ 250. 

Of course, in order to produce more accurate predictions for 
the observability, it is necessary to take into account the frequency 
of these events at high-z and the exact emission spectrum of the 
source: we defer this task to future work. 


4 DISCUSSION AND CONCLUSIONS 


In this work we have investigated the radiation-hydrodynamic evo¬ 
lution of the spherical accretion flow onto a DCBH with initial 
mass Mo = 10^ M© and gas mass in the parent halo Mgas = 
9.6 X 10® M©. 

The IMBH accretes for ~ 142 Myr with an average duty- 
cycle ~ 0.48 and on average super-critical accretion rates, with 
fsdd — 1.35, i.e. M, ~ 0.1 M© yr“^. The emitted luminosity in¬ 
creases with time, as a consequence of the progressive rise in the 
mass inflow. The radiation pressure creates strong waves moving, 
with velocities as high as ~ 20 kms“^, in the outer (r > 0.5 
section of the inflow. These waves produce shocks in some regions 
of the flow. 

At the end of the main evolutionary phase ~ 90% of the gas 
mass has been accreted onto the compact object, in spite of an av¬ 
erage super-Eddington emission, while ~ 10% has been ejected 
with outflows. The accretion is terminated when the emitted lumi¬ 
nosity reaches the value ~ 3 x 10^®ergs“^ ~ SLsdd and the 
related radiation pressure ejects all the remaining gas mass (which, 
at the flnal time, is a factor of ~ 25 lower than the initial one). 
We estimate that this flnal burst of radiation, lasting at least 5 yr in 
the rest-frame, should be observable by the future JWST. We have 
identifled three different phases of the accretion (the initial phase, 
the main accretion phase and the flnal burst), detailing their main 
characteristics in turn. 

We predict that the accretion flow, on average, occurs at mildly 
super-critical rates for the total evolution of the system (except 
the very initial transient phase). Recently, [Alexander & Natarajan| 
< |2014| ), [Volonteri & Silk| ( [20T4] ) and |Madau et al.| ( |2014| l, but see 


Volonteri & Rees 


2005| ), have suggested that brief but strongly 


also 

super-critical accretion episodes (with rates as large as fsdd > 50) 
might explain the rapid black hole mass build-up at 20 > z > 7. 
Very large and prolonged accretion rates may be sustainable in the 
so-called “slim disk” solu tion ([Begelman & Meier|1982[|Paczyn-| 


ski & Abramowicz||1982[[Mineshige et al.|2000[|Sadowski||2009 

201 1| |, an energy-advecting, optically thick flow that generalizes 
the standard thin disk model ( |Shakura & Sunyaev|1976| l. In these 
radiatively inefficient models the radiation is trapped (see also the 
diffusion condition in Appendix) into the high-density gas and is 
advected inward by the accretion flow: this happens when the pho¬ 
ton diffusion time exceeds the time scale for accretion. This would 
allow a very mild dependence of the emitted luminosity L on the 
normalized accretion rate fEdd, which is usually described as a log¬ 
arithmic relation ( [Mineshige et al.||2000[ |Wang & Netzer|[2003] ): 
LjLEdd ~ 2[1 + ln(/£;dd/50)], valid for fEdd > 50. 

In our work, we analyze the accretion flow on very large scales 
(R Rb, i.e. the accretion disk is beyond our resolution limit) and 
photons are never trapped (see the Appendix). For this reason, the 
accretion flow is radiatively efficient (LjLEdd oc fEdd) and the 


^ Estimate performed with the JWST prototype Exposure Time Calculator 
(ETC): http://jwstetc.stsci.edu/etc 


accretion rates are only mildly super-critical: fEdd — 1.35 on av¬ 
erage. In our case, the idle phases are caused by the necessity to 
re-establish the downward accretion flow after the radiation burst, 
while in the strongly super-critical models they are caused by the 
need for replenishing the gas reservoirs (e.g. by galaxy mergers, 
see e.g. | Volonteri & Silk||2014| ). In our simulation, at smaller ra¬ 
dial distances (R <C i^s) we expect that the radiation eventually 
reaches the trapping condition. This is neglected in the present im¬ 
plementation of the simulation, but could critically modify the ra¬ 
diative properties of the source, especially in the light of the afore¬ 
mentioned recent studies. Some aspects of the simulation may be 
improved, namely: 

(i) The accreting gas, at smaller radii, should form an accretion 
disk. If the accretion flow has a non-zero angular momentum with 
respect to the central body, the gas will reach a centrifugal barrier 
(caused by the steeper radial scaling of the centrifugal acceleration, 
~ with respect to the gravitational acceleration, ~ r“^) from 
which it can accrete further inward only if its angular momentum is 
transported away. This would at least partly modify the irradiation 
mechanism. 

(ii) A full spectral analysis of the source needs a more accurate 
description of the interaction between radiation and matter. 

(hi) At smaller radial distances the photons should be trapped 
and the accretion flow should become energy-advective, i.e. radia¬ 
tively inefficient, as described above. An appropriate modeling of 
the inefficient accretion flow would then be required. 

(iv) The magnetic held may also signiflcantly affect the accre¬ 
tion flow structure and behavior, as already pointed out in e.g. |Sad-| 
|owski et ^ ( |2014| ) and [McKinney et aL] ( |2014| ). The inclusion of 
an appropriate modeling of the magnetized plasma would then be 
required. 

This work is the basis upon which a full spectral analysis of 
these sources is to be constructed: this would be the key to unveil 
the eventual existence of IMBHs during the Cosmic Dawn era. The 
existence of IMBHs at high redshifts, although not yet confirmed 
by observations, would represent a breakthrough in our knowledge 
of the primordial Universe. Their presence would be also relevant 
for at least two reasons. First, the formation of IMBHs in the early 
Universe would ease the problem of the presence of SM BHs with 
masses M, ~ 10^ M© at redshifts as high sls z ^ 7.085 { [Mortlock 
jet al.|20lT] ). These massive seeds would play a role of paramount 
importance in giving a jump start to the accretion process. Second, 
the formation of IMBHs at high redshifts could provide a possi¬ 
ble interpretation of the near-infrared cosmic background fluctua¬ 
tions ( |Yue et alT]|2013| ) and its recently detected cross-correlation 
with the X-ray background ( jCappelluti et al.|2013| l. This interpre¬ 
tation would be even more plausible if the primordial population of 
IMBHs is proved to be highly obscured. The observation of IMBHs 
could then provide the missing pieces for the solution of these in¬ 
triguing puzzles. 


We thank L. Ciotti, G. S. Novak and M. Volonteri for useful 
comments and suggestions. 


5 APPENDIX 

This Appendix contains some more technical details about the mod¬ 
ules HD and RT, such as the boundary conditions for velocity 
and density, the two-stream approximation, the heating and cool¬ 
ing terms and the photon diffusion condition. 
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5.1 Additional Boundary Conditions 

In addition to the boundary conditions for the luminosity, already 
detailed in Sec. we need to specify the behavior of velocity and 
density at the innermost and outermost cells of the spatial grid. The 
spatial boundary conditions are: (i) outflow for the inner boundary 
(with restrictions on velocities, see below) and (ii) void for the outer 
boundary. An outflow boundary condition forces the derivatives of 
the quantities of interest to be zero, i.e. artiflcially extends the spa¬ 
tial domain outside the boundary. The restriction for the boundary 
velocity Vb is the following: 

I ij"min ) if U ( 7 * 777 ,177,) 0 

Vb — \ (jo) 

if7;(rmzn)>0 

and is meant to prevent the replenishment of the computational do¬ 
main by the gas coming from unresolved spatial scales. A void 
boundary condition, on the contrary, constrains the quantity of in¬ 
terest to be zero outside the computational domain: the system com¬ 
posed by the IMBH and its parent halo (up to ~ 2Rb) is isolated 
in space. 


5.2 The two-stream approximation 

The RT method we used is based on [Novak et al.| ( [20T^ and relies 
on the two-stream approximation, i.e. the luminosity is expressed 
as the sum of an ingoing and an outgoing radiation stream. When 
the optical depth is low, photons of the ingoing stream at any given 
radius ro are likely to successfully traverse the inner parts of the 
halo (r < ro) and emerge as outgoing photons at the same radius, 
but with (f = (/? + TT if is the azimuthal angle. In this case, all 
the radiation emitted by a source term is to be added to the outgo¬ 
ing stream. If, instead, the optical depth is large, the ingoing pho¬ 
tons are likely to be absorbed for r < ro. Then, only half of the 
emitted photons should be added to the outgoing stream. The other 
half should be added to the ingoing stream, where they will in due 
course be absorbed. The resulting equations for the two radiation 
streams are: 

= 47rrVr - pKLout (37) 


dLin 

dr 


47rr^(l - '0)^ - pnLin 


(38) 


where ^(r) is the fraction of photons emitted at a given radius that 
are likely to belong to the outgoing stream. The simplest estimate 


1 

1 


rf 

2 

_1 + e-^_ 


max(rf, r^) 


(39) 


where r is the optical depth from r to inflnity and ri is the radial 
distance from the center where r reaches unity. These equations 
must be complemented with the boundary conditions: 


Lin {Vmax) 
Lout (vmin^ 



(40) 


and the expression for the radiative acceleration, Eq. must be 
changed into the following one: 


(^rad 


L^i^Lout Lin) 


(41) 


5.3 Heating and Cooling 


This paragraph deals w ith the no n-adiab atic regime, i.e. the 


source/sink term E (Eqs. 25|26 

and 37|38 1 . Differentlv from the 

implementation in Novak et al. 

2012 | 

I, where E accounts for the 


energy transfer among different frequency bands, in our case the 
interaction between matter and radiation is purely elastic (i.e. the 
frequency of the interacting photon is unchanged) and this term ex¬ 
presses the energy emitted or absorbed by matter per unit time and 
per unit volume. 

For a gas at T < 10^ K (i.e. below the atomic hydrogen cool¬ 
ing threshold) the energy equation is purely adiabatic. Instead, for a 
gas at T > 10^ K, the term E takes into account the bremsstrahlung 
cooling and the atomic cooling (recombination and collisional ex¬ 
citation, for H and He): 

S = 2 - 0 XiounS (42) 


Here, n is the number density and xjj is the fraction of photons emit¬ 
ted at a given radius that are likely to belong to the outgoing stream 
(see above in this Appendix). The term 2 (tjj — |) is the fraction 
of transmitted (i.e. not absorbed) photons, since ^ is deflned as 
'0 = (1 + ptrans/ 2 ), whcrc ptrans is thc transmission probabil¬ 
ity. The term S includes all the coefficients for brem sstrahlung and 
atomic cooling (in units ergcm^ as reported in Maselli et al. 
|200^ and |CiK[T^ . 


Moreover, the term Xion is the ionized fraction, which 
accounts for the fraction of atoms that can contribute to 
the bremsstrahlung. This last term is calculated as Xion = 
{ab/x + 1 )^ where ab is the recombination rate and 7 is the col¬ 
lisional ionization rate. Both of them are expressed in cm^s“^ and 
are deflned as (see e.g. [Maselli et al.|2003| ): 


ab ■■ 


8.4 X 10“ 


T " 

1000 J 


1.0 + 


71V 

106 ) 


(43) 


7=1.27xl0-''v7e-i”'/^ 




(44) 


5.4 Photon Diffusion Condition 

Given the very high values of the optical depth (as high as ~ 100 
for the FS simulation), it is important to check whether a proper 
treatment for photon diffusion is required. Photons in the diffusive 
regime are advected inward with the gas, rather than being diffused 
out of the accretion flow. In this condition, the trapped infalling 
material should be considered a “quasi-star”, similar to the phase 
advocated by [Begelman et al.[p008[ t, with an atmosphere in Local 
Thermodynamic Equilibrium: the emission from the inner section 
of the accreting flow would then be thermal. 

[Begelman| ( [1978[ ) gives a very practical way to assess the oc¬ 
currence of the diffusive behavior of photons. Photons displaced at 
some radius r are trapped-in if: 

T(r) > 1 (45) 

c 

Throughout the paper, we refer to this formula as the diffusion con¬ 
dition, which is never met in our grid, because the maximum values 
reached are of order r{r)v(r)/c ^ 10“^. This proves that the pho¬ 
tons inside the simulated spherical volume in our work are never in 
the diffusive regime. It is likely that this condition is actually met at 
smaller radii, where the gas may form a structure quite similar to a 
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Stellar atmosphere. We defer the investigation of this possibility to 
future work. 
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